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1  Objectives 

The  objective  of  our  effort  was 

•  to  develop  mathematical  algorithms  and  corresponding  efficient  nu¬ 
merical  simulation  tools  for  modeling  propagation  of  acoustic  waves 
through  the  human  head, 

and,  subsequently, 

•  to  apply  these  tools  to  simulate  energy  transfer  to  the  inner  ear,  and 
to  assess  the  noise  induced  damage  to  the  human  hearing  system. 

2  Overview  of  the  approach 

It  is  a  well  established  fact  that  integral-equation  formulations  provide  the 
most  accurate  solutions  to  wave  problems.  They  require,  however,  solv¬ 
ing  dense  systems  of  linear  equations.  Traditional  methods  of  solving  such 
systems  (e.g.,  through  matrix  inversion)  are  characterized  by  the  computa¬ 
tional  complexity  and  memory  requirements  of  the  order  of  AT3.  Here  N 
is  the  number  of  unknowns,  or  a  number  of  volume  elements,  necessary  to 
resolve  acoustic  wave  variations  as  well  as  geometrical  details  of  the  ob¬ 
ject.  For  complex,  anatomically  realistic  models,  N  can  easily  become  of 
the  order  of  tens  of  millions.  Hence,  despite  their  reliability  and  accuracy, 
the  conventional  integral-equation  based  methods  become  computationally 
prohibitively  intensive  to  provide  solutions  for  realistic  problems  of  interest. 

Dining  the  last  fifteen  years  a  significant  progress  has  been  made  in  the 
development  of  fast  frequency-  and  time-domain  integral-equation  solvers 
and,  as  a  result,  the  ability  of  accurate  and  fast  numerical  simulations 
of  wave  propagation  and  scattering  in  complex  media  has  dramatically 
improved.  Matrix  compression  techniques,  such  as  the  Adaptive  Integral 
Method  (AIM)  l,  2],  the  Fast  Multiple  Method  (FMM)[3],  the  pre-corrected 
FFT  method[4],  the  sparse-matrix  canonical  grid  method  (SMCG)[5],  have 
been  developed.  They  allow  solving  large  linear  sets  of  equations  with  dense 
matrices  utilizing  storage  and  execution  times  characteristic  of  problems  in¬ 
volving  sparse  linear  systems.  The  physical  idea  behind  such  methods  is  that 
evaluation  of  interactions  at  large  distances  requires  less  resolution  than  at 
small  distances.  Consequently,  the  computational  complexity  and  memory 
requirements  of  the  compression  methods  scale  approximately  linearly  with 
the  number  of  unknowns  N . 
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The  underlying  element  of  our  approach  is  the  Fast  Fourier  Transform 
(FFT)  based  AIM  matrix  compression  method,  initially  developed  in  the 
context  of  electromagnetics  for  solving  large-scale  problems,  and  described 
in  detail  in  Ref  [2].  Adaptation  of  this  formulation  to  large-scale  acoustic 
problems  was  the  initial  step  of  our  effort. 

The  main  reason  for  choosing  the  FFT-based  compression  method  rather 
than  other  compression  techniques  is  that  it  provides  superior  efficiency  in 
the  treatment  of  both  volumetric  problems  and  sub-wavelength  (tetrahe¬ 
dron  size  much  smaller  than  the  wavelength)  discretizations.  We  note  that 
sub-wavelength  geometry  regions  constitute  dominant  portion  of  the  head 
geometry  model. 

The  final  goal  envisioned  by  us  when  undertaking  this  approach  was 

•  to  develop  an  efficient  simulation  tool  capable  to  solve  large-scale  prob¬ 
lems  without  compromising  accuracy  of  the  solution  (due  to  the  use 
of  an  integral  equation  based  solver  with  non-lossy  impedance  matrix 
compression), 

•  to  develop  capability  to  work  with  geometries  characterized  by  realistic 
geometrical  details  and  material  properties, 

and  thus  to  provide  a  numerical  tool  for  modeling  multiple  physical  mech¬ 
anisms  of  acoustic  energy  transfer  to  the  human  ear  and  of  assessing  the 
noise  induced  damage. 

3  Summary  of  the  results 

Our  work  comprised  the  following  main  efforts: 

•  development  of  fast  integral  equation  based  solver  for  acoustic  wave 
propagation  through  inhomogeneous  media, 

•  development  of  fast  volumetric  integral  equation  formulation  for  solv¬ 
ing  high-contrast  acoustic  problems  (e.g.,  biological  tissues  embedded 
in  air), 

•  a  parallel,  distributed  memory  implementation  of  the  fast  acoustic 
integral  equation  based  solver, 

•  extension  of  the  developed  approach  to  elasticity. 
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The  results  of  our  work  have  been  presented  in  two  articles  published  in 
the  Journal  of  Acoustical  Society  of  America  [8,  7].  Two  additional  papers, 
A  cousto- elastic  integral  equation  based  numerical  simulation  tools  for  anal¬ 
ysis  of  sound  wave  interactions  with  human  hearing  system  and  design  of 
high-noise  protection  devices  and  A  parallel  distributed-memory  implemen¬ 
tation  of  a  fast  acoustic  integral- equation  solver  are  in  the  process  of  being 
submitted  for  publication. 

The  Sections  below  give  the  brief  summary  of  the  work  performed. 

3.1  Fast  volumetric  integral  solver  for  acoustic  wave  propa¬ 
gation  through  inhomogeneous  media 

We  adaptated  the  Fast  Fourier  Transform  (FFT)  based  AIM  matrix  com¬ 
pression  method, [1,  2]  initially  developed  in  the  context  of  electromagnetics, 
to  solving  large-scale  volumetric  integral  equations  problems  in  acoustics. 

In  the  AIM  approach  the  original  basis  functions  are  expanded  into  sets 
of  auxiliary,  fax-field  equivalent,  point  sources  located  on  nodes  of  a  regular 
Cartesian  grid.  As  the  result,  the  far-held  part  of  the  interaction  matrix 
becomes  a  Toeplitz  matrix  (or,  more  precisely,  a  product  of  a  Toeplitz  and 
sparse  transformation  matrices).  The  Toeplitz  property  allows  computation 
of  the  matrix- vector  product  by  means  of  FFTs  with  the  N  log  N  solution 
complexity.  The  near-held  part  of  the  interaction  matrix  is,  by  construction, 
sparse.  Hence,  its  contribution  to  the  solution  complexity  is  of  the  order  of 
O(N). 

In  comparison  to  other  methods,  we  hnd  the  AIM  fast  solution  scheme 
particularly  suitable  for  solving  acoustics  problems.  Its  main  advantage  is 
that  it  can  be  applied  to  a  wide  frequency  range  -  from  tens  of  Hz  to  tens 
of  kHz  -  with  only  a  minor  modification  of  the  matrix  compression  scheme 
from  the  low-  to  high-frequency  mode.  (In  the  FMM  approach,  on  the 
other  hand,  entirely  different  algorithms  have  to  be  used  for  low  and  high 
frequencies.  The  p re-corrected  FFT  method  is  less  efficient  than  AIM  at 
higher  frequencies.  Finally,  for  the  considered  applications,  the  cost  of  the 
SMCG  method  is  significantly  higher  due  to  the  necessity  of  computing  a 
number  of  derivatives  of  the  Green  function.) 

The  developed  algorithm  makes  possible  simulations  involving  realistic 
geometries  described  in  terms  of  a  few  million  unknowns,  characterized  by 
highly  sub-wavelength  details,  and  large  density  contrasts. 

Results  of  our  work  are  presented  in  detail  in  [8].  Examples  involving 
calculations  of  acoustic  field  distribution  in  a  human  head  described  in  terms 
of  approximately  a  million  unknowns  are  also  shown  there. 
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3.2  Two-stage  acoustic  integral-equation  solver 

Conventional  Lippmann-Schwinger  integral  equations,  when  applied  to  high- 
contrast  problems  (e.g.  problems  involving  acoustically  dense  objects  (bi¬ 
ological  tissues)  immersed  in  a  lowr-density  medium  (air)),  may  potentially 
become  ill-conditioned. 

An  important  feature  of  such  high-contrast  problems  is  the  fact  that 
most  of  the  incident  wave  energy  is  reflected  from  interfaces  between  the 
low-  and  the  high-density  media,  due  to  a  large  impedance  mismatch  be¬ 
tween  the  materials.  This  circumstance  does  not  cause  computational  prob¬ 
lems  in  solving  purely  scattering  problems,  when  these  are  formulated  in 
terms  of  the  boundary-element  method,  i.e.,  as  surface  integral  equations 
involving  boundary  conditions  on  the  surface  of  the  scatterer.  Difficulties, 
however,  arise  if  one  attempts  to  determine  pressure  distribution  within 
an  inhomogeneous  body  by  means  of  the  volumetric  integral  equations  (the 
Lippmann-Schwinger  equations).  In  this  case,  the  fact  that  only  a  small 
fraction  of  the  incident  energy  penetrates  inside  the  scatterer  causes  the 
solution  inside  the  object  to  be  poorly  determined.  Mathematically,  the  in¬ 
tegral  equations  describing  the  system  become  ill-conditioned:  the  solution 
is  then  strongly  dependent  on  small  changes  in  the  incident  field,  and  may 
become  completely  unreliable. 

We  developed  a  solution  method  which  addresses  the  encountered  diffi¬ 
culties.  In  the  proposed  approach  an  ill-conditioned  high-contrast  problem 
is  solved  in  two  stages:  (1)  in  the  first  stage  we  solve  the  surface  integral 
equation  (as  in  boundary-element  methods,  although  using  volumetric  ele¬ 
ments);  (2)  in  the  second  stage  ,  we  solve  the  volume  problem  with  a  modified 
incident  field,  defined  in  terms  of  the  original  field  and  the  solution  of  the 
surface  problem.  Both  the  surface-  and  the  modified  volume  problems  are 
well  conditioned.  The  procedure  is  rigorous,  it  does  not  involve  expansions 
in  the  ratios  of  the  material  parameters,  and  it  does  not  require  alternat¬ 
ing  solving  the  surface  and  the  volume  equations  (although  each  of  the  two 
problems  may  be,  and  usually  is  -  for  large  systems  -  solved  iteratively). 

The  approach  is  described  in  detail  in  [7]. 

3.3  Distributed-memory  parallelization  of  the  acoustic  integral- 
equation  solver 

As  we  stated  above,  due  to  the  utilization  of  the  FTT  based  compression 
scheme,  the  computational  complexity  and  memory  requirements  of  our 
solver  scale  approximately  linearly  with  the  number  of  unknowns  N.  This 
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allows  for  routine  solutions  of  problems  up  to  about  1,000,000  unknowns  on 
a  single  processor  system  equipped  with  4  GB  memory. 

However,  in  order  to  perform  realistic  large-scale  simulations,  involving 
tens  of  millions  of  unknowns,  development  of  a  parallel  distributed-memory 
(MPI)  code  version  which  could  run  on  a  PC  cluster  with  typical  storage  of 
about  2  GB  memory  per  processor  becomes  essential. 

Distributed  memory  parallelization  of  an  integral  equation  based  solver  is 
highly  nontrivial,  since  the  nature  of  the  algorithm  gives  rise  a  large  amount 
of  inter-processor  communication.  The  three  main  computational  stages  to 
be  considered  in  the  parallelization  process  are: 

(i)  geometry  processing  and  distribution  of  the  geometry  data, 

(ii)  construction  of  the  stiffness  matrix, 

(iii)  the  iterative  solution  scheme. 

A  scalable  parallel  implementation  of  the  solver  requires  that  practically 
all  geometry  data  are  distributed  across  processors  with  a  minimal,  if  any, 
replication.  At  the  same  time,  during  the  construction  of  the  stiffness  ma¬ 
trix,  processors  must  have  access  to  some  global  geometry  data  in  order  to 
evaluate  near-field  couplings  between  sources  and  field  variables  assigned  to 
different  processors.  We  solve  this  problem  by  partitioning  the  geometry 
into  non-overlapping  slices  and  by  assigning  them  to  different  processors. 
However,  in  the  matrix  construction  stage,  we  temporarily  assign  several 
adjacent  slices  (of  combined  thickness  equal  at  least  to  the  near-field  range) 
to  each  processor  and  “weld”  them  into  a  single  stack.  This  welding  pro¬ 
cedure  allows  us  to  treat  the  entire  stack  of  slices  as  a  complete  geometry, 
without  having  to  introduce  any  additional  connectivity  data  relating  geom¬ 
etry  elements  in  adjacent  slices  (which  would  have  been  necessary  if  slices 
were  treated  as  separate  geometries). 

The  most  essential  part  of  the  algorithm  is  the  FFT-based  stiffness  ma¬ 
trix  compression  and  the  associated  fast  matrix-vector  multiplication  pro¬ 
cedure.  The  fact  that  Fourier  transforms  have  to  be  evaluated  in  all  three 
spatial  directions  requires  a  global  rearrangement  (“transposition”)  of  the 
data  and  hence  a  large  amount  of  inter-processor  communication. 

For  these  reasons  we  build  the  parallelization  scheme  of  the  code  around 
the  parallel  FFT  algorithm.  We  take  here  advantage  of  the  availability 
of  a  widely  used  FFTW  package  [9],  which  allows  operations  on  FFT  data 
distributed  spatially  (in  the  form  of  “slices”)  across  the  processors.  Its  MPI- 
parallelized  implementation  is  available  in  both  the  current  version  3  and 
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the  previous  version  2.  We  opted  for  the  more  recent  MPI  alpha  version 
3.2alpha3  (recently,  in  November  2008,  this  version  has  been  replaced  by 
3.3alphal). 

However,  a  literal  application  of  the  FFTW  routines  would  have  been,  in 
our  case,  inefficient.  The  equivalent  cartesian-grid  representation  of  sources 
and  fields,  used  in  our  fast  compression  algorithm,  requires  zero-padding 
(with  the  index  range  extension  by  the  factor  of  two  in  each  direction)  in 
order  to  eliminate  aliasing  in  the  discrete  Fourier  transformations,  the  con¬ 
ventional  FFTW  routines  would  locally  operate  on,  and  globally  exchange, 
zero-filled  buffers.  Our  application  of  the  FFTW  package  allows  a  balanced 
distribution  of  the  padding  storage  across  the  processors  and  reduces  the 
amount  of  communication  (the  latter  by  avoiding  sending  and  receiving  zero 
or  irrelevant  padding  data).  Such  an  efficiency  could  not  have  been  achieved 
by  applying  a  single  global  FFT  transformation  from  the  package. 

In  spite  of  requiring  a  large  amount  of  communication,  the  parallel  FFT 
implementation  in  our  solver  achieves  a  nearly  perfect  scaling  with  the  num¬ 
ber  of  processors  P  (typically,  a  speedup  by  a  factor  about  P2/3  instead  of 
the  ideal  speedup  ~  P).  We  tested  this  behavior  for  up  to  several  hundred 
processors. 

The  parallelization  procedure  implemented  in  the  present  code  assumes 
a  tetrahedral  geometry  discretization.  However,  it  is  applicable,  with  only 
minor  modifications,  to  other  types  of  geometry  discretization  (e.g.,  to  node- 
based  basis  functions). 

In  Section  4  below  we  present  some  large  scale  simulations  obtained  with 
the  parallel  solver.  We  show  solutions  for  a  human  head  model,  with  and 
without  a  helmet,  and  with  various  materials  used  as  padding  in  the  space 
between  the  head  and  the  interior  helmet  surface.  The  geometry  models 
were  discretized  with  approximately  10  million  tetrahodra. 

The  parallelization  procedure  is  described  in  detail  in  the  attached  draft 
of  a  manuscript,  and  will  be  submitted  for  publication. 

3.4  Construction  of  the  elasto-acoustic-solver 

Our  subsequent  research  effort  focused  on  the  construction  of  a  full  elasto- 
acoustic  integral  equation  formulation  and  on  the  development  of  the  related 
solver.  In  our  work,  we  built  on  our  experience  acquired  in  the  construction 
of  a  purely  acoustic  solver:  in  particular  when  overcoming  difficulties  asso¬ 
ciated  with  large  contrast  cases  and  strongly  sub-wavelength  contributions. 

Three  candidate  versions  of  the  solver  were  constructed: 
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-  first  order  coupled  acousto-elastic  solver  for  the  4-dimensional  un¬ 
known  vector  composed  of  three  components  of  the  displacement  field 
u(r)  and  the  pressure  field  p(r), 

-  second  order  acousto-elastic  solver  for  the  3-dimensional  unknown  vec¬ 
tor  composed  of  three  components  of  the  displacement  field  u(r), 

-  first  order  coupled  acousto-elastic  solver  for  the  9-dimensional  un¬ 
known  vector  composed  of  three  components  of  the  displacement  field 
u(r)  and  six  (out  of  total  nine)  independent  components  of  the  stress 
tensor  field  f(r). 

The  above  formulations  differ  in  the  relative  content  of  either  first  or  sec¬ 
ond  order  derivatives  of  material  parameters,  and  in  the  treatment  of  high 
contrast  contributions.  At  this  point  the  choice  of  the  numerically  optimal 
formulation  poses  still  an  open  question  and  more  realistic  calculations  are 
needed  to  be  carried  out  to  answer  it. 

In  addition  to  the  listed  above  three  candidate  versions  of  integral  equa¬ 
tions  formulations  we  also  constructed  a  first  order  coupled  acousto-elastic 
solver  applicable  to  geometries  composed  of  piecewise  homogeneous  regions. 
Such  version  of  the  solver  involves  exclusively  unknown  values  of  the  traction 
field  on  interfaces  between  regions  characterized  by  different  Lame  param¬ 
eters  A,//  and  the  material  density  p.  The  main  reason  for  devising  such  a 
method  and  constructing  a  surface  integral  equation  solver  was: 

-  to  be  able  to  model  highly  complex  elements  of  the  inner  ear  with 
surface  elements, 

-  to  quickly  assess  the  importance  of  mechanisms  contributing  to  energy 
transfer  to  the  cochlea  within  the  framework  of  a  simplified  but  rel¬ 
atively  realistic  geometry  model  we  constructed  on  a  parallel  STTR 
effort  (we  briefly  describe  the  model  in  Appendix  A  below), 

-  to  verify  the  validity  of  the  purely  volumetric  integral  equation  solver 
(in  addition  to  the  analytic  solution  for  a  multi-layered  elastic  sphere 
we  constructed  on  a  parallel  STTR  effort). 

As  we  noted,  our  choice  of  the  optimal  integral  equation  representation 
to  be  implemented  in  the  final  version  of  the  elasto-acoustic  solver  remains 
still  an  open  issue.  However,  we  derived  numerically  expressions  for  matrix 
elements  appearing  in  all  of  the  candidate  versions  of  the  elasto-acoustic  inte¬ 
gral  equations.  In  our  derivations  we  used  a  new  method  of  evaluation  of  in¬ 
tegrals  of  kernel  elements  involving  dipole  terms  of  the  Green  function.  The 
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procedure  offers  both  analytical  simplicity  and  numerical  accuracy.  It  does 
not  require  singularity  extraction  procedure  since  it  reduces  six-dimensional 
volumetric  integrals  to  four-dimensional  surface  integrals  with  nonsingular 
integrands. 

The  derived  expressions  are  applicable  to  linear  basis  functions  on  tetra¬ 
hedral  and  triangular  supports. 

The  article  version  of  the  approach  described  above  is  in  preparation. 


4  Examples 

Pressure  distributions  in  amultilayer  sphere  -  comparison  with  the 
analytical  solution.  As  the  first  example  we  present  comparison  of  the 
results  obtained  with  our  solver  with  analytical  solutions.  Figs.  1  and  2 
show  results  for  the  pressure  distribution  in  a  layered  sphere,  with  the  outer 
two  shells  chosen  to  represent  skin  and  bone,  and  the  interior  of  the  sphere 
described  by  mechanical  parameters  of  the  brain  tissue.  The  outer  radius 
of  the  sphere  is  14  cm,  the  thickness  of  the  skin  layer  is  1  cm,  and  so  is  the 
thickness  of  the  bone  layer.  The  calculations  employed  a  tetrahedral  mesh 
with  the  tetrahedron  sizes  of  about  2.8  mm;  hence  the  thickness  of  the  shell 
was  about  three  times  the  tetrahedron  size.  The  total  number  of  tetrahedra 
is  about  N  =  5.5  million. 

Fig.  2  shows  distribution  of  the  absolute  value  of  the  pressure,  |p(r)|,  for 
the  incident  wave  of  frequency  /  =  3  kHz). 


Figure  1:  Discretization  of  the  layered  sphere  with  N  =  5.5  million  tetrahe¬ 
dra. 
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Figure  2:  Distribution  of  the  absolute  value  of  pressure,  p(r)|,  on  the  plane 
passing  through  the  center  of  the  layered  sphere,  computed  analytically  (left) 
and  by  means  of  the  described  acoustic  solver  (right),  for  the  discretization 
with  N  =  5.5  million  tetrahedra. 

The  high  quality  agreement  between  the  analytical  results  and  code  pre¬ 
dictions  for  this  high-contrast  problem  provides  a  support  for  the  accuracy 
of  our  code. 

Pressure  distributions  in  a  human  head  model  in  the  presence  and 
absence  of  a  helmet.  The  next  example  involves  a  realistically  shaped 
model  of  a  human  head  and  a  model  of  helmet. 

We  assumed  a  steel  helmet  and  the  space  between  the  helmet  and  the 
head  filled  either  with  air  or  with  cork.  In  this  computation  we  assumed  the 
head  model  filled  with  a  homogeneous  material,  with  mechanical  properties 
of  bone.  The  models  are  placed  in  the  field  of  a  harmonic  acoustic  wave 
of  frequency  5  kHz,  incident  laterally  on  the  right  ear  of  the  head  model. 
The  corresponding  wavelength  in  air  is  A  =  6.8c.m.  Acoustic  parameters  of 
materials  used  in  this  and  similar  computations  are  listed  in  Table  1. 

Geometries  were  discretized  with  tetrahedron  sizes  (edge  lengths)  of 
about  3  mm,  resulting  in  (a)  N  ~  2,700,000  tetrahedra  for  the  head,  and 
(b)  iV  ~  4, 700, 000  for  the  head  and  helmet  system.  The  computations  were 
carried  out  on  a  Linux  cluster  with  the  InfiniBand  interconnect  network,  on 
108  processors  for  the  model  (a)  and  128  processors  for  the  model  (b).  The 
total  computation  times  were  (a)  about  50  minutes  and  (b)  about  2  hours, 
the  longer  time  for  the  case  (b)  due  mostly  to  the  larger  number  of  iterations 
in  the  solution.  One  iteration  required  4.6  s  in  the  problem  (a)  and  6.6  s  in 
the  problem  (b).  We  estimate  that,  in  both  cases,  the  overall  computation 
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time  can  be  reduced  by  30  to  50  %  by  optimizing  the  matrix  construction 
stage  of  the  code. 


Table  1:  Acoustic  properties  of  materials  used  in  the  simulations 


material 

p/pO 

rr 

bone 

1777.0 

0.1524 

brain 

835.4 

0.0564 

cork 

150.0 

0.65 

steel 

6667.0 

0.0035 

Figure  3:  Pressure  distribution  on  the  surface  of  the  head  in  the  absence 
of  the  helmet  (left),  pressure  distribution  on  the  surface  of  the  helmet  and 
the  head  with  the  space  between  the  head  and  the  helmet  filled  with  air 
(center),  pressure  distribution  on  the  surface  of  the  head  and  the  helmet 
with  the  space  between  the  head  and  the  helmet  filled  with  cork  (right). 

Fig.  3  shows  pressure  distribution  on  the  surface  of  the  head  in  the 
absence  of  the  helmet,  as  well  as  pressure  distributions  on  the  surface  of  the 
helmet  and  the  head  with  the  space  between  the  head  and  the  helmet  filled 
with  air  and  cork  respectively.  In  all  cases  we  observe  a  surface  wave  being 
formed  on  the  geometry  exterior. 

In  Fig.  4  we  show  pressure  distribution  inside  the  human  head  model  in 
the  presence  and  in  the  absence  of  a  helmet. 

The  results  show  a  nontrivial  behavior  of  the  solutions  and  exhibit  phys¬ 
ical  phenomena  which  may  be  relevant  in  the  design  of  protective  devices. 

In  the  case  of  the  head,  Fig.  4(a),  the  pressure  is  maximal  at  the  entrance 
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to  the  ear  canal,  and  it  is  relatively  smoothly  distributed  inside  the  head. 
In  fact,  the  solution  is  suggestive  of  a  resonance-type  (P-wave)  behavior: 
the  pressure  changes  sign  along  the  approximately  vertical  line  seen  in  the 
Figure. 

The  solution  for  head  and  helmet  system,  Fig.  4(b),  is  quite  different. 
It  exhibits  a  distinct  oscillatory  behavior  along  the  surface  of  the  helmet 
and  in  the  region  filled  by  cork.  This  region  appears  to  have  properties  of  a 
“waveguide1.  Because  of  the  cork  density  being  significantly  lower  than  that 
of  the  surrounding  materials  (the  helmet  and  the  head),  and  the  resulting 
impedance  mismatch  at  the  boundaries,  the  wave  tends  to  be  trapped  in 
that  region.  Since  the  refractive  index  of  cork  is  not  much  different  from 
that  of  air,  wave  oscillations  are  relatively  rapid.  We  stress,  however,  that 
the  physical  picture  suggested  by  Fig.  4(b)  would  change  if  we  included 
dissipative  (attenuation)  effects  in  the  filling  material,  e.g.,  if  we  considered 
a  strongly  damping  porous  material  characterized  by  a  complex  refractive 
index. 

We  also  note  that,  for  the  particular  frequency  considered  here,  the  pres¬ 
ence  of  a  helmet  completely  changes  the  pressure  distribution  in  the  head, 
but  does  not  reduce  its  maximum  value  (the  data  in  Fig.  4(a)  and  Fig.  4(b) 
are  plotted  in  different  scales). 


Figure  4:  Pressure  distributions  in  the  coronal  plane  for  (a)  the  human  head 
model  and  (b)  the  system  consisting  of  the  human  head  and  a  steel  helmet 
models,  with  the  in-between  space  filled  by  cork.  The  models  are  subject  to 
an  acoustic  wave  of  unit  pressure  amplitude  and  frequency  5  kHz,  incident 
from  the  left.  The  maximum  pressure  value  is  about  4  in  (a)  and  15  in  (b). 
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5  Appendix  A:  Test  model  of  a  human  head 

On  a  parallel  STTR,  Phase  II  effort,,  we  constructed  a  test  model  of  a  human 
head.  The  model  consists  of  a  skin,  skull  and  brain  tissue,  with  the  cochlea 
embedded  in  the  skull.  In  addition,  we  also  constructed  a  model  of  a  helmet 
and  a  padding  material  placed  between  the  head  and  the  interior  of  the 
helmet  surface.  The  components  of  the  head  model  are  shown  in  Figs.  5. 
The  helmet  model  is  shown  in  Fig.  6. 

The  above  elements  constitute  a  minimal  yet  relevant  set  of  geometry 
components  required  in  carrying  out  simulations  of  energy  deposited  within 
the  cochlea  region  and,  subsequently,  in  the  determination  of  auditory  effects 
of  the  propagating  wave.  The  helmet  geometry  allows  us  to  investigate  the 
influence  of  the  padding  material,  as  well  as  of  possible  air  gaps  between  the 
head  and  the  helmet  structure  on  the  energy  distribution  inside  the  head 
model. 


Figure  5:  The  external  ear,  the  cochlea,  and  the  skull  models  used  in  the 
simulations. 
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Figure  6:  Model  of  the  external  head  surface  and  the  helmet.  The  arrow 
indicates  the  direction  of  the  incident  wave. 
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